Effects of nanoparticles and surfactant on droplets in shear flow 
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We present three-dimensional numerical simulations, employing the well-established lattice Boltz- 
mann method, and investigate similarities and differences between surfactants and nanoparticles as 
additives at a fluid-fluid interface. We report on their respective effects on the surface tension of 
such an interface. Next, we subject a fluid droplet to shear and explore the deformation properties 
of the droplet, its inclination angle relative to the shear flow, the dynamics of the particles at the 
interface, and the possibility of breakup. Particles are seen not to affect the surface tension of the 
interface, although they do change the overall interfacial free energy. The particles do not remain 
homogeneously distributed over the interface, but form clusters in preferred regions that are stable 
for as long as the shear is applied. However, although the overall structure remains stable, individual 
nanoparticles roam the droplet interface, with a frequency of revolution that is highest in the middle 
of the droplet interface, normal to the shear flow, and increases with capillary number. We recover 
Taylor's law for small deformation of droplets when surfactant or particles are added to the droplet 
interface. The effect of surfactant is captured in the capillary number, but the inertia of adsorbed 
massive particles increases deformation at higher capillary number and eventually leads to easier 
breakup of the droplet. 

PACS numbers: 47.11.-j 47.55. Kf, 77.84.Nh, 
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I. INTRODUCTION 

Stabilizing emulsions by employing nanoparticles is a 
very attractive tool in the food, cosmetics, oil and med- 
ical industries. This method of emulsification comple- 
ments the traditional use of surfactants - amphiphilic 
molecules - as emulsification agents. Using nanoparti- 
cles can have many advantages, such as reduced cost and 
toxicity and the possibility of tailor-made nanoparticles, 
which may include useful properties other than being an 
emulsifier, such as ferromagnetic particles [1] or Janus 
particles [2| . Although the effects of both emulsifiers can 
be similar, the underlying physics is very different d, Q. 

Amphiphiles are chemical compounds which have both 
hydrophilic and hydrophobic properties, restricted to 
specific groups of the molecules. For example, surfac- 
tants are characterized by their hydrophilic "head" and 
hydrophobic "tail(s)". When they are located at the 
fluid-fluid interface the possibility exists for both parts of 
the molecule to reside in their preferred fluid. This makes 
it energetically favourable for them to accumulate at the 
interface, with a distinct alignment. This process lowers 
interfacial tension and prevents the demixing of two im- 
miscible fluids. As such, it gives rise to the possibility 
of complicated structures, such as micelles and lamel- 
lae, gyroid mesophases and the aforementioned emulsion 
droplets [BHH. 
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Nanoparticles also find it energetically favourable to 
adsorp to a fluid- fluid interface, however, this happens for 
a different reason. Maintaining such an interface requires 
more energy per unit area than maintaining a particle- 
fluid interface, and the adsorption of a particle removes 
the former. Because of the scale of the energy differ- 
ences involved - orders of magnitude larger than thermal 
fluctuations - this adsorption process tends to be irre- 
versible [3[. In this way, neutrally wetting particles do 
not affect surface tensions directly, but only change the 
interfacial free energy. 

When such particles are used to stabilize an emulsion of 
discrete droplets of one fluid suspended in another, con- 
tinuous, fluid, the result is known as a "Pickering emul- 
sion 

|a,II3|. The particles in these mixtures block Ost- 
wald ripening, which is one of the main processes leading 
to drop coarsening in emulsions. Hence, blocking this 
process allows for a long-term stabilization of such an 
emulsion. They are also a source of complex rheology 
due to the irreversible adsorption of the particles as well 
as interface bridging because of particle monolayers [TTJ— 
fl3j . More recently, the use of nanoparticles has led to the 
discovery of the "bicontinuous interfacially jammed emul- 
sion gel" (commonly referred to as "bijel" ), first predicted 
by numerical simulations [3] and later confirmed exper- 
imentally 0, EH. In a bijel, an interface between two 
continuous fluids (as opposed to having separate droplets 
of one fluid) is covered and stabilized by particles. The 
effect of parameters such as fluidifluid ratio and particle 
wettability on the final phase a demixing system trans- 
forms into has been investigated numerically [l7H20j. 

The differences between the behaviour of amphiphiles 
and nanoparticles and between their underlying mechan- 
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ics as described above ensure that many properties of 
systems including nanoparticles cannot be explained by 
theories based solely on the physics of amphiphiles. For 
nanoparticle-stabilized systems, new models have been 
developed (and verified experimentally), which take into 
account the features of these systems that have no ana- 
logue in surfactant systems, such as the contact angle 
of the nanoparticles, strong capillary forces between the 
particles or the pH value and electrolyte concentration of 
the solvents [21[. Quantitatively, however, the descrip- 
tion of these systems still leaves to be desired. 

To properly understand the behaviour of large-scale 
mixtures with many complex interfaces, such as Picker- 
ing emulsions and bijels, one first needs a fundamental 
understanding of the processes involved on smaller scales. 
Research was performed to understand in detail how the 
presence of a nanoparticle 1221 or the collective behaviour 
of multiple nanoparticles [23, 24] affects a flat interface. 
In this work we investigate the stabilizing effect of am- 
phiphiles or hard spherical nanoparticles on curved inter- 
faces, modeled by a single droplet of a fluid suspended in 
another fluid. 

Droplets subjected to shear flow display many kinds 
of interesting behaviour, such as deforming away from 
a spherical shape, exhibiting an inclination angle with 
respect to the shear direction and breaking up into 
smaller droplets (beyond a critical capillary number) [25| . 
Nanoparticles adsorped at the droplet interface show 
an inhomogeneous distribution and a non-trivial motion 
over the droplet surface. Their presence also affects the 
deformation and inclination properties of the droplet. We 
study all these effects in detail in the current article. 

Computer simulations are a valuable tool to compare 
these systems directly, and we choose to employ the lat- 
tice Boltzmann (LB) method, which is well-established in 
the literature (cf. [26]), for our research. The LB method 
is an alternative to traditional Navier- Stokes solvers, and 
extensions have been developed to allow for multiple flu- 
ids and their interactions [2^-32] , amphiphiles 0, [33| 
and finite-sized particles of arbitrary shape and wetta- 
bility which can interact with the fluids as well as each 
other ESS. 



In section |TT] we introduce the simulation method in 
detail. Section IIIII reports and explains our findings on 
surface tensions in systems of a droplet stabilized by sur- 
factant and nanoparticles. The behaviour of the parti- 
cles adsorped to the droplet interface when the droplet is 
subjected to shear is discussed in section HVl This is fol- 
lowed by a an analysis of the effect of nanoparticles and 
surfactant on the deformation properties and inclination 
angles of these droplets. The breakup of droplets is then 
briefly discussed. Finally, conclusions and an outlook are 
provided in section [Vj 



II. SIMULATION METHOD 
A. The lattice Boltzmann method 

The lattice Boltzmann method has proven itself to be 
a very successful tool for modeling fluids in science and 
engineering (gfj [39|, H(| . Compared to traditional Navier- 
Stokes solvers, the method allows an easy implementa- 
tion of complex boundary conditions and - due to the 
high degree of locality of the algorithm - is well suited 
for implementation on parallel supercomputers @, [Io| . 

The method is based on the Boltzmann equation, with 
its positions x discretized in space on a cubic lattice with 
lattice constant Ax and with its time t discretized with 
a timestep At: 

/?(x + c t At, t + At) = /?(x, t) + fi?(x, t), (1) 

where /f(x, t) is the single-particle distribution function 
for fluid component c, being propagated over the lattice 
with a discrete set of lattice velocities and 

a c (v A _ /P(x,t)-/T q (p C (x > t),U C (x > t)) 

is the Bhatnagar-Gross-Krook (BGK) collision opera- 
tor [4l|. Here, / 2 eq (p c ,u c ) is the third-order equilibrium 
distribution function 
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t c is the relaxation time for component c and Q are the 
coefficients resulting from the velocity space discretiza- 
tion [42]. We use a three-dimensional lattice and a 
D3Q19 implementation (i = 1, . . . , 19), which is to say 
that Axi = CiAt connect a lattice site with its near- 
est neighbours and next-nearest neighbours on the lat- 
tice. The Navier- Stokes equations can be recovered from 
Eq. (pQ). The macroscopic densities are given by p c (x, t) = 
p c (x,t)/pQ = /?(x, t), with /?q being a reference den- 
sity for component c. For clarity of notation, the tilde 
is omitted from the densities from now on. The macro- 
scopic velocities are u c (x, t) = ff (x, t)ci/p c (x, t) in 
the low Knudsen number and low Mach number limit. 
The speed of sound on the lattice is 
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from which one can calculate the kinematic viscosity of 
a fluid component as 



(5) 



dAti At-2 



For convenience, the lattice and time constants are taken 
to be Ax = At = 1 from now on. In all simulations 
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presented here, we have chosen r c = 1 for all components, 
which then implies v c = 1/6 for all components. The size 
of the simulation volume is denoted as V^o* = n x -n y -n z . 



B. Multicomponent lattice Boltzmann 

When further fluid species d with a single-particle dis- 
tribution function ff (x, t) are to be modeled, an interac- 
tion force Fq(x,£) is calculated locally according to the 
approach of Shan and Chen [27j : 

F^(x,t) = -* c (x,t)^ 5cc ^vl/ c '(x',t)(x'-x) , (6) 



with g cc f a coupling constant and \I/ C (x, t) a monotonous 
weight function representing an effective mass. Through- 
out this work, this function takes the form 



M/ c (x,t) = vl/(p c (x,t)) = l-e 



-P°(x,i) 



(7) 



This force is then incorporated into the collision term 
in Eq. (pQ) by adding to the velocity u c (x, t) in the 
equilibrium distribution the shift 



Au c (x,t) 



r c F^(x,t) 
P c (x,t) 



(8) 



Furthermore, this forcing affects the macroscopic bulk 
velocity as 

u c ( t) = Es/f(M)c, _ \ ¥ c ( t y (9) 

v ' ; P c (x,t) 2 cv ' ) w 

In our case, the coupling strength g cc > is negative in or- 
der to obtain de-mixing and the sum over x' in Eq. (|6j) 
runs over all sites separated from x by one of the discrete 
velocities In the binary fluid systems we refer to the 
fluid of the droplet (d) and the medium (m) as "red" 
fluid (r) and "blue" fluid (6), respectively. To simplify 
statements about the fluid:fluid ratio on lattice sites, we 
introduce the order parameter 0(x, t) = p r (x, t)— p 6 (x, i), 
referred to as "colour". The LB method is a diffuse in- 
terface method, with an interface width of typically 5 
lattice sites, depending weakly on the coupling strength 
gbr- Owing to this, there will typically also be a small but 
non-zero density of red fluid population in the medium 
and of blue fluid population inside the droplet. This will 
be touched upon in greater detail in section IIII Al 



C. Amphiphiles 

Amphiphiles can be introduced to LB simulations in 
various ways. While Benzi et al. have proposed a method 
that can reach even vanishingly low surface tensions by 
including mid-range interaction forces [43|, Hif , we avoid 
taking into account additional Brillouin zones and in- 
stead follow a model proposed by Chen et al. 0, l33l EBf. 



Although this method suffices to recover the qualitative 
behaviour of surfactant, it is limited in the surface ten- 
sion reduction it can effect - 60% reduction being the 
largest achieved in our simulations. However, availabil- 
ity of larger reduction was deemed unnecessary for the 
purpose of the present work. 

In addition to having its own set of distribution func- 
tions //(x, t), the amphiphilic surfactant (s) has a dipole 
vector d(x, t) associated with it, representing the aver- 
age orientation of the amphiphiles at a lattice site. The 
direction of this dipole vector can vary continuously. Its 
propagation is given by 

/ s (x,i+l)d(x,i+l) = 

£(//(x-Ci,*)d(x-c;,i)). (10) 

i 

Here, the tildes denote the post-collision values - for a 
quantity Q\\ Q\ = Q\ -f . The relaxation of the dipole 
vector can also be described by a (vector) BGK process 
as 

w x w \ d(x,t) -d eq (x,f) , . 

d(x, t) = d(x, t) - V ; ^ V (11) 

with r d the relaxation time of the dipole orientation to- 
wards a local equilibrium d eq (x, t). Furthermore, the 
force terms as described in Eq. (j6j) are extended to ac- 
count for the forces the amphiphiles exert on the red and 
blue fluids: 



F c (x,t)=F£(x,t)+F§(x,t) , 



(12) 



where the lower indices denote the source of the force and 
C and S refer to "colour" and "surfactant" , respectively. 
The new addition to the force term takes the form 

F c s (x, t) = -2^ c (x, t)g C8 d(x + c„ t) • 0^ s (x + c, , t) , 

(13) 

where g cs is the force coupling constant between an ordi- 
nary and the amphiphilic species and Oi is a second-rank 
tensor defined as 



Oi = 1-3- 



(14) 



with 1 the second-rank identity tensor. Similarly, the 
forces acting on the amphiphiles can be split into contri- 
butions from amphiphiles and ordinary fluid: 



F s (x,t) = F£(x,t) + F^(x,t). 
These take the forms 



(15) 



F£(x, t) = 2* s (x, t)d(x, t) -J29csJ2 ^* C ( X + c *^) 

(16) 



and 

F|(x,t) = -^2 fe ^ s (x,t)- 

ll c ill 

^r(x + c„t)(d(x + c„t)^ r d(x,t)c, 

i ^ 

+ [d(x + a, t)d(x, t) + d(x, i)d(x + a,t) 



(17) 



respectively. The coupling constant g ss should be neg- 
ative to model attraction between two amphiphile tails 
and repulsion between a head and a tail. For a full deriva- 
tion of these equations, cf. jfj. 



D. Nanoparticles 

Nanoparticles are discretized on the lattice and cou- 
pled to both fluid species by means of a modified bounce- 
back boundary condition as pioneered by Ladd [33 - 
13(1 [46[ , resulting in a modified lattice Boltzmann equa- 
tion 

/?(x + Ci,t + l) = /f(x + Ci ,i) + fif(x + Ci ,*) + C , (18) 

where C is a linear function of the local velocity of the 
particle surface, and i are defined such that = — cj. 
Wherever x is occupied by a particle, Eq. (pQ) is replaced 
by Eq. ([T8|) . The particle configuration is evolved in 
time, solving Newton's equations in the spirit of classical 
molecular dynamics simulations. As the total momen- 
tum has to be conserved, an additional force acting on 
the particle is needed to compensate for the momentum 
change of the fluid caused by Eq. (ITS]) : 



F(t) = (2/f(x + c i) t) + C)ci. 



(19) 



As the simulation evolves in time and a particle moves 
around, the configuration of lattice sites occupied by the 
particle changes. When a site is newly occupied by a 
particle, the fluids on that site are deleted and their mo- 
mentum is transferred to the particle through a force 



F(f) 



p c (x,t)u c (x,f). 



(20) 



Lattice sites which have been newly vacated by a particle 
also have to be treated. In Ladd's original algorithm for 
a single fluid, the initial fluid density pf nit would be used. 
However, in the case of a multicomponent system this 
would cause artefacts, in particular for the case of parti- 
cles adsorped to an interface: fluid b would be initialized 
where only fluid r ought to be present and vice versa. To 
prevent such problems from occuring, a density 
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FN 



V(x + C ipN ,£), 



(21) 



is defined, averaged over the Afn neighbouring fluid lat- 
tice nodes x iFN , separated from x by the velocity vector 



c^ FN . The fluid on the vacated site is initialized with 
populations 

/?(x,t) = C(x,t) • /; q (u surface (x,t),p new (x,t)), (22) 

where u sur f ace (x, t) is the local velocity of the particle 
surface. Due to non-zero repulsive Shan-Chen forces be- 
tween the particle surface and the surrounding fluid, the 
effective fluid density close to the particle surface might 
be slightly smaller than the bulk density leading to a 
mass drift over time if one chooses p new (x, t) = p c (x, t). 
To suppress this effect we apply a correction which keeps 
the total mass constant on long time scales, with small 
fluctuations (of the order of 10 -4 of the total mass) on 
shorter time scales [18j : 



Pnew( X ^) = p C (*,t) (l - Co 



Pfnit ^>ox 



(23) 



where Ap c (t) is the total mass error of color c at time 
t, and Co can be used to tune the strength of the cor- 
rections. In this work, Co = 2500 is used. To prevent 
instabilities, we restrict this density to be not larger or 
smaller than the highest and lowest surrounding density, 
respectively. 

The potential between the particles is a Hertz potential 
which approximates a hard core potential and has the 
following form for two spheres with identical radii r p j!?} : 



(24) 



Yj 1 1 is the 



§h = K H (2r p -rijY for 

and zero otherwise. Here, = ||r^ 
distance between the centres for two spheres z, j located 
at and r^, respectively, and Kh is the force constant, 
which we choose to be Kh = 100. Apart from the di- 
rect interaction described by the Hertz potential we cor- 
rect for the limited description of hydrodynamics when 
two particles come very close by means of a lubrication 
correction. If the number of lattice points between two 
particles is sufficient - at least one fluid site - the LB 
algorithm reproduces the correct lubrication force auto- 
matically. If particles approach beyond this limit, the 
flow is no longer sufficiently resolved. The error can be 
corrected by an additional force term 



-rplub 
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jTij ■ (Ui-Uj) 



1 



2r n 



1 

a: 



(25) 

with Ui and Uj the velocities of particles i and j, respec- 
tively and Tij the unit vector pointing from the centre of 
particle i to the centre of particle j [36|. Furthermore, 
we choose a cut-off of this lubrication force A c = 2/3. 

The force in Eq. (j6j) also includes interactions between 
a lattice node outside of a particle and a lattice node 
inside a particle. To calculate these interactions the lat- 
tice nodes in the outer shell of the particle are filled with 
a "virtual" fluid corresponding to the density defined in 
Eq. (f2T|) : p^ irt (x, t) = p c (x, t). This density is assigned 
to the population density / r c est (x, t) for which c rest = 0. 
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Advection and collision are not applied to this virtual 
fluid. 

A system of two immiscible fluids and particles is con- 
sidered. We define a parameter Ap, the particle colour, 
which allows to control the interaction between the parti- 
cle surface and the two fluids. If Ap has a positive value, 
we add it to the red fluid component as 

P ; irt = f + Ap. (26) 

Otherwise we add its absolute value to the blue fluid as 



^virt 



Ap. 



(27) 



By changing Ap it is possible to control the contact angle 
P of the particle. The dependence of the contact angle 
on the particle colour can be fitted by a linear relation, 
where the slope depends on the actual simulation param- 
eters. A particle colour Ap = corresponds to a contact 
angle of P = 90°, i.e. a neutrally wetting particle. For 
a more detailed description of our simulation algorithm 
the reader is referred to [l8[ . 

An alternative method to introduce solid particles to 
free-energy-based multicomponent LB simulations was 
introduced by Stratford et al. [3, [Hj|, while Joshi and 
Sun published applications of the multiphase Shan- Chen 
model with suspended particles [37j. 



E. Boundary conditions 

The simulation volume is bounded at the x = 1 and 
x = n x planes by Lees-Edwards shear boundary condi- 
tions [49| , which avoid spatial inhomogeneities that occur 
when shear is induced by moving walls. These boundary 
conditions have been adapted for use in LB simulations 
by Wagner and Pagonabarraga [Eoj], and the reader is 
referred to this publication for technical details. In our 
simulations the boundary conditions are set up in such a 
way as to effect a shear rate 7 = 2u s / (n x — 1) in the z- 
direction. The remaining sides of the system are subject 
to ordinary periodic boundary conditions 0, H[ • 

III. SURFACE TENSION 
A. Theory 

The Young-Laplace equation relates the pressure dif- 
ference AP over the interface between two fluids to the 
surface tension a: AP = — crV-n, with n the surface nor- 
mal. For a spherical (undisturbed) droplet of one fluid 
of radius Rd inside another fluid this equation takes the 
form 



a = 



R d AP 



(28) 



Calculating the correct pressure jump AP = P^ — 
P m > over the interface, where P^ is the pressure in- 
side the droplet and P m is the pressure in the medium, 



requires some care. For a single-component and single- 
phase system, local pressure in LB can be calculated us- 
ing the simple relation P(x) = c|p(x) (here and in all fu- 
ture equations, the time dependence has been suppressed 
in our notation). However, when using the multicompo- 
nent Shan-Chen model for a ternary system - consisting 
of simple fluid species red r and blue b and a surfactant 
species s - there is a non-zero presence of the local mi- 
nority fluid throughout the system and we have to use 
the more complicated expression 

^W(x)+p fc (x) + /9 s (x)+ 

+ ^ 5cc ^ c (x)* c '(x)+ 5ss * s (x)* s (x), (29) 

c/c' 

which takes into account pressure contributions of the 
fluid- fluid interactions. 

Because of the diffuse interface in LB simulations one 
has to make sure that the measurements are performed 
far enough away from the interface, so that the density is 
(almost) constant in the neighbourhood. We have veri- 
fied that the density profiles of the system in equilibrium 
are flat on the inside and outside of the droplet as little 
as five lattice sites away from the isosurface where the 
colour field is zero. Hence, this effect does not cause a 
problem in these cases. We therefore take a spatial av- 
erage of the pressure in the centre of the droplet over a 
small neighbourhood (5 3 cube of lattice sites) as P^, and 
the local spatial average in a corner of the system (which 
due to the periodic boundary conditions is the furthest 
away one can get) as P m . Densities of the fluids c - de- 
noted as p c d and p^ - can now be defined in a similar 
manner. 

Calculating the radius of the droplet is also non-trivial, 
again due to the diffuse interface. We have investigated 
three distinct approaches, whose results have been in 
agreement up to less than a lattice site - two methods 
based on detection of the <fi = isosurface and one based 
on total mass and density of the red fluid. The latter 
method has been chosen, since it can most easily be ex- 
tended to the case of added particles, which will be ex- 
plained below. We consider the idea that all the surplus 
population of the red fluid ought to be contained in a 
sphere of constant density. We define a local effective 
density p r eE (x) = p r (x) — p^ to account for the non-zero 
density of red fluid outside of the droplet. This effective 
density is used to calculate the total droplet mass 



(30) 



See Fig. [T] for an illustration of this process. 

Using the relation for the droplet volume Va = 
Md/{p r d — p r m ) and assuming sphericity of the droplet 
leads to 



M d 



47r / Pd-Pm- 



(31) 
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FIG. 1. A ID profile of the local densities p r (x) (solid curve) 
and p b (x) (dashed curve) along the z-axis and centered in the 
x — y plane, as used in the calculation of the droplet mass Ma- 
The droplet density of the red fluid p r d and medium densities 

are taken at the center (circle) and edge (squares) of the 
domain, respectively. The shaded area denotes the summed 
effective total density of red fluid, which takes into account the 
non-zero "background" density of the red fluid in the medium 

by subtracting it from the local densities. 



When nanoparticles are adsorped at the droplet inter- 
face (which could change the shape of the = isosur- 
face dramatically depending on the number of particles 
and their position, validating the choice of this particular 
method), a correction term is needed to account for these 
particles. Recalling that the radius of the spherical parti- 
cles is denoted r p , we define a new effective volume of the 
droplet VJ ff = Vd + V p , and approximate V p ~ ^ (if r p)> 
where n p is the number of particles, expressing that we 
expect half of the particle volume to be on the inside of 
the interface of the droplet, adding its volume to the vol- 
ume derived from the red fluid. Thus, the final equation 
for the radius of the droplet is given by 



Rd 



3 

47T 



M d 



P r 3 



(32) 



From Eq. ([29]) and Eq. ([28]) one can see that the mea- 
sured surface tension depends on the fluid densities - 
linearly in first order, but in a more complicated fashion 
in the cross terms, where the form of the effective mass 
function ^ plays a role (cf. Eq. (|ZJ)). In light of this, we 
keep the initial density of the simple fluid species con- 
stant across simulations. 



B. Effect of amphiphiles 

We now proceed to study the effect of added surfac- 
tant on the system. The system is initialized as follows: 



FIG. 2. Main plot: plotting the pressure jump over the fluid- 
fluid interface against the inverse droplet radius, the Young- 
Laplace equation for a spherical droplet allows to calculate 
surface tensions as the slope of fitted linear functions going 
through the origin AP = a (2/Rd) for various values of the 
Shan-Chen interaction parameter gb r - Inset: the surface ten- 
sion a is a monotonically increasing function of gbr- The 
lines connect the results taken from the linear fits, while the 
symbols show the averaged results of direct evaluation of the 
surface tension for single points on the curves. Direct calcula- 
tion of the surface tension is an accurate and efficient method 
when considering many systems with different parameters, 
which would require extra simulations with multiple droplet 
radii otherwise. 



a cubic simulation volume = 64 is consid- 

ered, and the initial droplet is chosen to have a radius 
of R 1 ^ = 0.3n x = 19.2 and is placed in the centre of 
the system. These values were chosen after determining 
the effect of the resolution of the lattice on the surface 
tension. The total variation when increasing the system 
size from 48 3 to 128 3 was seen to be less than 4%, and 
the best balance between accuracy and computational ef- 
fort was attained at 64 3 . This error is smaller than the 
errors from the sources described above. After discretiza- 
tion, the interior droplet sites are set to have a density 
P r = Pirnt an d P b — 0- Conversely, the medium sites 
have p r = and p b = pf nit , while the interface is crudely 
modeled by a linear density gradient over 5 lattice sites. 
Because of stability reasons and the shape of the effec- 
tive mass function ^ c we use pT nit = p^ nit = 0.7 in all 
results presented here. In the case of added surfactant, 
the density is set to p s = pf nit everywhere. The initial 
surfactant density varies from simulation to simulation 
and will always be reported explicitly. As the system ap- 
proaches its equilibrium state, surfactant accumulates at 
the interface, causing the local density at the interface to 
be higher by a factor of approximately two compared to 
the density in the bulk. Reaching the equilibrium state 
from this initialization can take a long time - to obtain 
stable results for the surface tension the simulations have 
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to run for tens of thousands of time steps for the systems 
described in this paragraph (and up to several hundred 
thousand timesteps in the case of a system with particles, 
as will be described in section UlI C|) . 

Firstly, we are interested in determining the sur- 
face tension as a function of the fluid-fluid interaction 
strength g^ r in the case of a binary fluid system. We 
fix the coupling constants related to the surfactant to 
limit the parameter space of interest and choose g rs = 
9bs = 9ss = —0.005. As discussed in section HTCl these 
have to be negative to properly model the behaviour of a 
surfactant. The actual values are chosen for their stabil- 
ity. There are also some restrictions on our choice of g^. 
For gb r < 0.10 the fluids become miscible when surfac- 
tant with the properties specified above is added, lead- 
ing to ill-defined interfaces and droplets. Furthermore, 
choosing g^ r > 0.15 leads to numerical instabilities (5lj . 
We therefore consider the values 0.10 < g^r < 0.15, re- 
stricting reachable surface tensions. Rewriting Eqn. [28] 
as AP = a (2/Rd) allows to extract a by considering it 
to be the slope of the pressure difference plotted against 
twice the inverse droplet radius. Linear fits through the 
origin correspond very well to the simulation results for 
0.10 < g br < 0.15 (cf. Fig. El- From this it follows that 
g\y r can be mapped onto the surface tension: a = cr(#& r ), 
with (J {gbr) a monotonically increasing function. The 
inset of Fig. [2] shows that calculating surface tensions di- 
rectly using a single droplet radius together with Eqn. [28] 
is an accurate and efficient method that does not require 
multiple simulations for a single choice of g^. 

The qualitative result of creating a ternary system by 
adding an amphiphilic surfactant component to the bi- 
nary droplet system is as expected: increasing surfactant 



density from 



0.0 to p? 



0.15 and p? u 



0.25 



lowers the surface tension by 30 to 50 percent (cf. the 
inset of Fig. [3]). As mentioned in section HTCl this rela- 
tively modest reduction is due to limitations of the sur- 
factant model used for these simulations. It is, however, 
sufficient to highlight the differences between the effect 
of amphiphiles and nanoparticles. To find a quantita- 
tive relation between surfactant concentration fraction 
^ s = Pinit/ (pfnit + ^fnit) > interaction strength and sur- 
face tension, it is useful to define 



rel 



(- 



1 9br, 



(33) 



where cr = a((j) s = 0). Plotting this quantity as a func- 
tion of s , the data points collapse onto a universal curve, 
as shown in Fig. [3] This illustrates the fact that the ef- 
fect of the surfactant scales with the interaction strength 
between the two no n- amphiphilic fluid species. We can 
use this data to obtain another mapping: a = <j(gi, r ,(j) s ) 
for fixed interaction strengths involving the surfactant 
species. These mappings will later be used in determin- 
ing capillary numbers for systems of a droplet subjected 
to shear. 
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FIG. 3. Main plot: by rescaling the effect of surfactant to 
Acr re i = (<j/o"o — 1) gbr the curves for different values of fluid- 
fluid interaction strength g br can be made to collapse, illus- 
trating the fact that the effect of added surfactant scales with 
gbr- The error bars of the data points are too small to be 
visible at this scale. Inset: surface tensions as a function of 
gbr for various concentrations of surfactant pf nit . The lines 
are not a fit, but included only to guide the eye. This shows 
qualitatively that the surfactant lowers the surface tension, as 
expected. Again, the error bars are too small to be visible in 
this plot. 



C. Effect of nanoparticles 

Next, the case of added (spherical and monodisperse) 
nanoparticles is considered. The fraction x of the droplet 
surface removed by the adsorped particles is a parameter 
of interest. The excluded surface area due to one neu- 
trally wetting particle is a spherical cap whose area is 

2irR d (R d 



given by = 2irRd {^Rd — \J R^ — r^j , from which fol- 
lows that the total coverage fraction of a spherical droplet 
is given by 



At 



Rd 



X : 



2R d 



(34) 



Since we use a diffuse interface method, any suspended 
particles have to be of sufficient size compared to the in- 
terface width to resolve their interfacial properties. In 
practice, this means a typical spherical particle needs to 
have a diameter of at least 10 LB length units, while a 
spherical droplet should be larger than the particles by 
an order of magnitude. Allowing then sufficient room 
for the deformation of the droplets to take place with- 
out undue finite size effects, these calculations remain 
computationally challenging, even for the case of a single 
droplet and a highly efficient massively parallel simula- 
tion environment. In order to be able to let the droplet 
deform sufficiently in later simulations we also elongate 
the system in the direction of the shear flow (z-direction): 
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FIG. 4. Surface tension change as a function of particle 
droplet surface coverage x (top x-axis, circles) and surfac- 
tant volume fraction (j) s (bottom x-axis, squares). Here, ao is 
the surface tension for the purely binary system (i.e. x — 
and (j) s = 0, respectively) with otherwise identical parame- 
ters. For all cases gbr — 0.10; for the system with surfactant 
g rs — gbs — gss — —0.005 and for the system with particles 
r p = 5.0, m p = 524, and P = 90°. Introducing 25% volume 
fraction of surfactant into the system lowers the surface ten- 
sion by almost 60%, while particles affect it only very weakly. 
The slight drop in measured surface tension for moderate val- 
ues of x is caused by errors introduced in the calculation of 
the droplet radius due to anisotropic particle distributions on 
account of spurious currents at the droplet interface. 



n x = n y = 256, n z = 512. The droplet is initialized as 
described above, with initial radius ^ nit = 0.3 -n x = 76.8 
and we choose g^ r = 0.10. The particles have a radius 
r p = 5.0 and are neutrally wetting (0 p = 90°). Further- 
more, they have a mass m p = 524, which corresponds to 
a density p p = 1 (taken with respect to the lattice). They 
are initialized on a spiral running over the surface of the 
initial droplet from the north to south pole, resulting in a 
very uniform initial distribution of particles at low com- 
putational cost [52j | . When the system is allowed to get 
into its equilibrium state, however, some pattern forma- 
tion of the particles on the interface occurs, due to the 
occurence of spurious currents near the interface (as also 
observed in similar modeling of liquid- vapour systems by 
Joshi and Sun [37j ). This effect is negligible when the 
system is not stationary: the currents are much smaller 
than the effect of applying shear to the system, or, for 
example, the effect of droplet movement in the forma- 
tion of a Pickering emulsion. In either case the particle 
ordering due to the spurious currents is destroyed. 

Adding particles with the aforementioned properties 
does not affect surface tension at all - the presence of par- 
ticles at the interface only changes interfacial free energy 
directly by taking away energetically expensive fluid- fluid 
interface and replacing it with cheap particle-fluid inter- 



faces. To clarify this, consider the free energy term F a re- 
lated to the surface tension of the interface of the droplet 
D, 



F n = 



adA, 



(35) 



3D 



which integrates the surface tension over the interface. 
For simplicity, the surface tension is taken to be con- 
stant over the interface. There are now two possibilities 
to reduce this energy. The first is to reduce the surface 
tension a, which is the effect of added surfactant. Be- 
cause a > and the integration only pertains to the 
fluid-fluid interface, the second possibility is to reduce 
the area of integration dD, which is effected by adsorped 
particles. The particles also add energy to the system by 
means of the interfacial energy between the particle and 
either fluid, but as this energy per unit surface area is 
much smaller than the fluid-fluid surface tension, the net 
effect is still a reduction of the free energy. 

A comparison of the addition of surfactant and 
nanoparticles to a binary system can be seen in Fig. |4j 
Due to the anisotropic distribution of the particles on the 
interface, errors are introduced in the calculation of the 
droplet radius for intermediate values of lowering the 
measured surface tension by up to 3%. At higher the 
anisotropy disappears, and with it the calculated change 
in surface tension, which returns to its original value for 
X ~ 0.5. In the system with surfactant an identical value 
of gb r = 0.10 is used. Unlike adding particles, adding sur- 
factant lowers the surface tension (a 60% drop in surface 
tension for pf nit =0.25). 



IV. DROPLET IN SHEAR FLOW 
A. Theory 

The system of a droplet of a fluid suspended in another 
fluid is subjected to simple shear flow, which causes the 
droplet to deform (cf. Fig. [5j). To analyze this process 
we first define a set of dimensionless variables. The di- 
mensionless deformation parameter 



D 



L-B 
L + B 



(36) 



introduced by Taylor [53|, |54| is used to describe the de- 
formation of the droplet, where L is the length and B 
is the breadth of the droplet. If the droplet is a perfect 
prolate ellipsoid the length and breadth can be related to 
the long and short axes, respectively, but in other cases 
a length and breadth of a more irregular shape can still 
be determined. One can easily see that for a spherical 
droplet L = B, hence D = 0, and for a strongly de- 
formed droplet, L ^> B, D — >• 1. Extraction of D from 
the data is effected through the symmetric moment of 
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FIG. 5. Representative deformation of a particle-covered 
droplet at Ca = 0.075, x — 0.55. The shaded planes at 
the top and bottom are subject to Lees-Edwards boundary 
conditions, inducing a shear rate 7 = 2u s / (n x — 1) in the 
z-direction, as discussed in section III El The shear causes 
droplet deformation D = (L — B) / \L + B) and an inclina- 
tion of the droplet of angle 9d, which is the angle the long 
axis of the droplet L forms with the shear direction z. 



inertia tensor 





"111 


h2 


hs 




I = 


I12 


I 22 


123 


(37) 
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I 23 


h3 





In order to define these moments of inertia, we first cal- 
culate the centre-of-mass position of the droplet 



X ' Pcom( X )> 



(38) 



where a cutoff density Pcutoff * s introduced to confine the 
summation to the droplet: 



with M el1 the mass of the ellipsoid and a, b and c the 
length of the axes. We now assume that the deformed 
droplet can be approximated by such an ellipsoid, and 
M el1 = M^ om . The set of equations obtained by combin- 
ing the eigenvalues of I with Eq. (|42j) can be solved for 
a, b and c. The length and breadth of the droplet are 
then defined as L = max(a, 6, c) and B = min(a, 6, c), 
respectively. 

A droplet thus deformed has lost its spherical shape 
and gains a preferred alignment. This is expressed 
through the inclination angle 64. It is calculated by tak- 
ing the eigenvector L corresponding to the long axis of 
the droplet of the moment of intertia tensor I, and cal- 
culating the arctangent of the quotient of its x and z 
components: 



arctan 



(43) 



A capillary number Ca can be defined as Ca = 
Hm'yRd/o', where /i m is the dynamic viscosity of the 
medium, 7 is the shear rate as imposed through the Lees- 
Edwards boundary conditions, Rd is the radius of the 
initial - undeformed, hence spherical - droplet and a is 
the surface tension. However, using this definition of the 
capillary number does not take into account the substan- 
tial distortion of the linear shear gradient caused by the 
presence of the droplet, which leads to a dependence on 
the size of the simulation volume, even in the case when 
only the resolution of the simulation is increased. A bet- 
ter characterization of the system can therefore be found 
in an effective capillary number: 



eff _ Mm7 eff #< 



Ca' 



(44) 



Pcom( X ) 



p r (x) if p r (x) > p^ utoff 
otherwise. 



(39) 



The cutoff density should fulfill the condition < 
^cutoff < Pd an d can be chosen freely within that range 
with negligible effect on the subsequent calculations. We 
use Pcutoff = 0-1 m this work. A droplet mass based on 
the density Pcom( x ) i s introduced for later use: 



M com 



^com( X )- 

xGVbox 



(40) 



Defining x = x — x^ om allows to express the elements of 
I as 



T ij = ^com( X ) (ll X H 2 Sij - XiXj ) , 

xeVbox 



where Sij is the Kronecker delta. The moment of inertia 
tensor I el1 of an ellipsoid of uniform density is a diagonal 
matrix with its non-zero elements given by 



ell 



M el1 



((l-(5 a )a 2 + (l-^ 2 )6 2 + (l-fe)c 2 ), 

(42) 



where an effective shear rate 7 eff is measured in the sim- 
ulation, instead of assuming the validity of an imposed 
shear rate set directly by an input parameter. Fig. [6] de- 
picts the measurement of 7 eff for a droplet with initial 
radius Rd = 39.2 in a system of size n x = n y = 128, 
n z = 256 and with u s = 0.05. Far away from the droplet 
7 eff & 7, but for those values of z over which the droplet 
extends, typically 7 eff > 7. The slope of the velocity 
gradient between the shear boundary and the droplet in- 
terface can be measured, which is then averaged over the 
length of the droplet to obtain the effective shear rate. 
This shear rate better characterizes the system. When 
the effective capillary number is used, taking into account 
the actual shear experienced by the droplet, the depen- 
dence of the deformation on the system size disappears, 
as shown in Fig. [7j where deformations of a droplet are 
plotted against both Ca and Ca eff . When the original 
capillary number is used, the deformation curves diverge 
as the system size increases from 64 2 • 128 to 128 2 • 256 
and 256 2 -512, while the curves collapse when plotted as 
a function of the effective capillary number. 

We also define the ratio of the droplet and medium vis- 
cosity A = i^d/l^m — 1 in all presented data, as well as a 
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FIG. 6. Representative z- velocity profiles of a droplet with 
initial radius Rd = 39.2, centered in a system with n x — 
riy = 128, n z = 256 and u s = 0.05. The cuts are taken in 
x-direction at y = 63 and taken through the droplet (z — 100) 
as well as at the edge of the periodic volume {z = 0). Also 
shown is what the imposed shear rate would look like in ab- 
sence of the droplet (7 is the slope of this line) . It is clear that 
far away from the droplet, the measured shear is almost undis- 
turbed and linear, while the droplet locally strongly disturbs 
the effective shear profile. We detect the droplet interface 
and calculate an effective shear rate 7 eff based on the slope 
of the profile in the region between the shear boundary and 
the droplet. Because of the deformation and inclination of 
the droplet this curve will generally not be symmetrical for 
the top and bottom shear boundaries for any particular given 
value of z, however, averaging over the length of the droplet 
restores this symmetry. 



Reynolds number Re ; 
number 

Re' 



and an effective Reynolds 



eff _ p m r s R 2 d 



(45) 



Due to the variation in system size and shear rate, the 
Reynolds number varies between approximately 0.6 < 
Re eff < 25, the effect of which will be discussed in sec- 
tion nvo 



B. Distribution of amphiphiles and nanoparticles 

To understand the effect of amphiphiles and nanopar- 
ticles on the deformation properties of the droplet, we 
first discuss how they position themselves at and move 
over the droplet interface as the droplet is sheared. 

The distribution of surfactant on a 2D cut through a 
sheared droplet is shown in Fig. [51 In this example, the 
shear rate is held constant at 7 = 0.002 and the initial 



FIG. 7. Dimensionless deformation D of a droplet in shear 
flow as a function of capillary number Ca = jjL T n A fRd/o' (left) 
and effective capillary number Ca eff = fimj eE Rd / & (right). 
Different symbols represent different system sizes. The cap- 
illary number computed from an assumed undisturbed shear 
profile gives rise to divergence in the relation between Ca 
and the deformation when the system size changes (lines are 
included to guide the eye). However, these points collapse 
on the curve which uses the effective capillary number, which 
takes into account the actual shear experienced by the droplet. 




FIG. 8. Distribution of surfactant in a system of size n x — 
riy = 64, n z = 128. The local surfactant density p s (x) is 
plotted on a 2D cut showing the centered x-z plane through 
a droplet sheared with constant velocity u s = 0.06 and 
P?nit — 0.25 (Ca eff = 0.16). The snapshot is zoomed into the 
droplet and as such does not accurately reflect confinement 
of the droplet or the elongation of the system. Surfactant 
accumulates at the droplet interface until saturation occurs. 
Compared to the remainder of the interface, a slightly higher 
local density is observed at the tips of the droplet (10 to 20 
percent). This is caused by convection of the surfactant [55[. 



the interface. When the system is subjected to shear, 
a slightly increased density of approximately 10 to 20 
percent is observed at the tips of the droplet, due to 
convection of the surfactant [55j . This behaviour is more 
0.25. As has been readily apparent for lower pf nit and is different from our 
mentioned in section UlI B( the surfactant accumulates at observations in the case of adsorped particles, as we will 



surfactant density is set to p? u 
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FIG. 9. Side-view examples of deformed droplets, for various 
particle coverage fractions: a) x — 0.00, b) % = 0.27 and c) 
X — 0.55. In these pictures the shear velocities are horizon- 
tal. In all these simulations gb r = 0.10, r p = 5.0, m p = 524, 
and P = 90°. One can see that although increasing x from 
to 0.27 does not strongly change the deformation of the 
droplet, the particles themselves do exhibit interesting be- 
haviour: they prefer to stay in the middle of the channel where 
the shear flow is weakest (recall that the top and bottom 
planes are moving inducing flow in opposite directions). This 
causes the formation of a band of particles near the equator 
of the droplet, with the axis through the poles in x-direction. 
For packings of higher density there is an interplay between 
shear forces and the curvature of the interface, which causes 
the aforementioned band to grow asymmetrically as the parti- 
cles prefer to occupy interface with high local curvature. The 
particles also exhibit tank-treading-like behaviour: they move 
around the interface following the shear flow. The combined 
effect of this tank-treading-like movement and the energy ar- 
guments described above lead to the formation of strings of 
single particles, being swept from the band near one tip to 
the other tip. 



show below. 

Even if the droplet interface is initially densely packed 
with particles, this will no longer be the case when the 
droplet deforms - the interfacial area increases while the 
number of particles remains constant. The particles then 
have freedom to move over the interface to some extent 
(cf. Fig. [9]). In all cases, however, detaching particles 
from the interface remains practically impossible. The 
particles are swept over the interface with increasing ve- 
locity as they move away from the centre plane of the sys- 
tem and up the shear gradient. If the particles would not 
be affected by the shear flow, they would prefer to occupy 
interface with high local curvature as can be explained 
by a geometrical argument: the interface removed by a 
spherical particle at a curved interface is larger than the 
circular area removed from a flat interface, and this ef- 
fect gets stronger as curvature increases. This explains 
why in this dynamic equilibrium, most particles can be 



FIG. 10. Main plot: normalized pair correlation function be- 
tween particles G n0 rm(r), for x — 0.41, r p — 5.0 and various 
capillary numbers. As the capillary number increases, the 
peaks both shift in position and increase in height. The for- 
mer effect is an indication of closer packing, while the latter 
corresponds to the observation of preferred regions for the 
particles as shown in Fig. [9]b). Inset: the height of the first 
three peaks of the normalized pair correlation function are 
shown as a function of the effective capillary number. The 
strongest effect is seen for the very first peak, which shows 
the largest growth in both relative and absolute sense and 
increases in height by almost a factor of 3. 



found at the tips of the droplet. This can be observed in 
Fig. [9] b) at high capillary number, where the relatively 
fiat sections of the interface at the top and bottom of the 
droplet have the lowest particle density and the strongly 
curved section of the interface near the centre plane is 
much more highly populated than the strongly curved 
section protruding farther into the shear flow. 

To quantify these phenomena, we employ a discrete 
pair correlation function 

G W = EE r^^llri-r.-H-^dfl, (46) 
i=l 3=1 Jr 

where S(x) is the Dirac delta function. Because of the 
system size, the domain of the pair correlation function is 
limited to < r < n x /2 = 128. Furthermore, we choose 
to employ G norm (r), which is proportional to G(r) and is 
normalized to tend to 1 as r tends to its maximum value. 
In Fig. [10] we show this normalized pair correlation func- 
tion for x = 0.41, r p = 5.0 and various capillary numbers. 
In the main plot, two features are readily apparent: the 
peaks of the function both shift in position and increase 
in height as the capillary number is increased. The for- 
mer effect indicates a denser overall packing of the par- 
ticles, which occurs despite the fact that more interfacial 
area becomes available as the droplet deforms. The latter 
effect corresponds to the emergence of preferred regions 
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FIG. 11. Rotational frequency of particles w as a function of 
their position on the y-axis for x — 0-55 and r p — 5.0 and 
various capillary numbers. The rotational frequency does not 
explicitly take into account the increase in the circumference 
of a cut of the droplet perpendicular to the ?/-axis due to de- 
formation. As the capillary number is increased through in- 
creased shear rate, the particles' rotation frequency increases, 
with particles in the middle of the droplet being both the 
fastest and getting the largest speedup, despite the fact that 
the particles in the middle travel the longest paths. The dif- 
ferences in frequency highlight that the movement is different 
from tank-treading as observed in, for example, vesicles [561 ]. 



FIG. 12. Main plot: deformation parameter D as a function 
of the effective capillary number Ca eff for various interaction 
strengths gb r and surfactant densities p s . The system size 
IS Tlx — 

= 64, n z — 128 and the surfactant interaction 
strenghs are fixed at g rs = gbs = g ss = —0.005. The capillary 
number is varied by changing the shear rate. Over the entire 
domain, the curves collapse onto a universal curve, and for 
small capillary number Taylor's law is recovered (dashed line). 
Note that the squares in this plot correspond to the squares in 
Fig. [71 Inset: The relation between the capillary and Reynolds 
numbers is a linear one. The slope is proportional to the 
surface tension (cf. Eq. (08])). 



for particles as described above. In the inset we show the 
height of the first three peaks of G norm (r). The first peak 
shows the largest increase, both in absolute and relative 
sense. This is caused by the fact that at x = 0-41 the 
band of particles around the droplet is not dense every- 
where, but is mostly restricted to patches near the tips 
of the droplet. Thus, particles having a close packing 
around them extending over more than one particle dis- 
tance (which would show peaks of higher order) are more 
rare than those with just closely packed neighbours. Fi- 
nally, we have observed that when this structure is es- 
tablished it is stable over time for as long as the system 
is subjected to a constant shear. When this shear is re- 
moved, the particles restore themselves to their former 
pattern, as described in section Ull C[ just as the droplet 
shape returns to that of a sphere. 

Even though the overall structure of the particles on 
the droplet interface remains stable over time, individ- 
ual particles move over the interface, performing a quasi- 
periodic motion. Their trajectories follow the motion of 
the shear flow and loop around the droplet with a rota- 
tional frequency uo. We demonstrate in Fig. [11] that this 
frequency is not constant for all particles, instead show- 
ing a dependence on the position of the particle along 
the y-cods. When deformation is considered for ellip- 
soidal cuts of the droplet along the ?/-axis, the deforma- 
tion is highest in the centre of the droplet, giving particles 



greater options for mobility that are also better-aligned 
with the shear flow, leading to increased particle veloci- 
ties. This is qualitatively different from the tank-treading 
behaviour observed in, for example, vesicles [56], which is 
characterized by a constant frequency for all points. We 
also observe that the average rotation frequency increases 
with increasing capillary number, in spite of the fact that 
the particles need to follow longer paths to complete one 
revolution as the droplet deforms. This increase in fre- 
quency is concentrated on the particles in the centre of 
the droplet, for the same reasons as mentioned above. 



C. Droplet deformation and inclination 

For small capillary number, Taylor predicts a linear 
dependence of the deformation of a droplet on the capil- 
lary number [HiJ Q , with a particularly simple form for 
equiviscous fluids (A = fid/l^m — 1) : 

^ = 19^ + 16^ Ca= g Ca 

16/i d + 16/i m 32 v } 

This law has been recovered in our simulations for the 
case of binary systems with various system sizes and 
interaction strengths, using the effective capillary num- 
ber introduced in section HV Al Combining Eq. ([44]) and 
Eq. (p[5]) one can derive a relation between the capillary 
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FIG. 13. Inclination angle 0d of a droplet as a function of 
the effective capillary number Ca eff for various interaction 
strengths g^r and surfactant densities p s . The system size 
IS Tlx — 

— 64, n z — 128 and the surfactant interaction 
strenghs are fixed at g rs = gbs = g ss = —0.005. The capil- 
lary number is varied by changing the shear rate. The high 
capillary numbers are only reached by using the surface ten- 
sion lowering effect of the added surfactant, hence the varying 
ranges of the datasets presented here. As the capillary num- 
ber increases, 0d can first increase beyond an angle of 45° 
(dashed line) due to viscous effects. Only when the forces 
induced by the shear start to dominate does the droplet align 
with the shear flow 15711. 



and Reynolds number: 



Re ett = a 



PmRd 
Mm 



Ca e 



(48) 



As we change the capillary number explicitly by changing 
the shear rate, the Reynolds number is proportional to 
the capillary number for a fixed value of the surface ten- 
sion. Inertial effects increase the deformation, thus the 
deformations at high capillary number are higher than 
predicted by the linear relation of Taylor. 

When a surfactant is added to the system, it lowers 
the surface tension of the interface, affecting the capillary 
number (but leaving the Reynolds number unchanged). 
Interaction strengths g^ r = 0.10 and g^r = 0.13 are used, 
while the surfactant interaction strengths are fixed at 
9rs — 9bs — Qss — 0.005 for the reasons mentioned in 
section UlI Bl The initial homogeneous surfactant densi- 
ties range from pf nit = 0.0 to pf nit = 0.3 in increments of 
0.05 and the system size is n x = n y = 64, n z = 128, with 
an initial droplet radius of ^ nit = 0.3 -n^ = 19.2. The de- 
formations for these systems are shown in Fig. [12j Since 
the change in surface tension directly enters the capil- 
lary number, all curves (including those not shown here 
for clarity) collapse onto a universal curve as a function 
of Ca eff , and Taylor's law is reproduced for small cap- 
illary numbers < Ca eff < 0.06. In the inset we show 



the relation between the capillary and Reynolds num- 
bers. It is clear that these relations are linear, the slopes 
are proportional to the surface tension (which is changed 
both implicitly and explicitly) , and agree with the values 
predicted by Eq. (|48j) . 

Inclination angles of the droplet in its steady state are 
plotted in Fig. [T3j In the case of Stokes flow, one would 
expect the inclination angle to be 45° for very small cap- 
illary number and to observe a decrease of this angle 
as the capillary number is increased, indicating a better 
alignment of the droplet with the shear flow. However, 
as inertia plays a role here we observe that in some cases 
0d first increases beyond 45°, before the inclination de- 
creases again and the droplet becomes elongated along 
the shear direction. When the steady inclinations are 
considered as a function of the Reynolds number, there 
exists a critical capillary number for which the inclination 
angle never exceeds 45°. Grouping the results at simi- 
lar capillary numbers as datasets, we estimate this to be 
Ca®f it w 0.11. Our observations are consistent with re- 
sults obtained by Singh and Sarkar, using a front-tracking 
finite- difference method [571 ]. 

We now consider a system with nanoparticles as ad- 
ditives. The fluid- fluid interaction strength is held fixed 
at gb r = 0.10. As before, the particles have a radius of 
r p = 5.0 and are neutrally wetting (0 P = 90°). Initially, 
we choose their mass to be m p = 524, as in section UlI CI 
As discussed previously, the introduction of finite-sized 
particles introduces a lower bound on how small the 
simulation volume can be to accomodate enough par- 
ticles on the interface and to avoid finite-size effects. 
For this reason, the simulation volume is chosen to be 
n x = n y = 256, n z = 512, with an initial droplet radius 
of ^ nit = 0.3 • n x = 76.8, still keeping it as small as 
possible to avoid excessive calculation time. The number 
of particles is varied as n p = 0, 128, 256, 320, 384, 446 
and 512, which results in a surface coverage fraction of 
X = up to x — 0.55. Again, the capillary number is 
changed by changing the shear rate. Some examples of 
the deformations thus realised are shown in Fig. [9j for 
Ca eff = 0.04, 0.08, 0.12 and X = 0.0 (a), X = 0.27 (b) 
and x — 0-55 (c). 

Although the effect of addition of surfactant on the de- 
formation and inclination of the droplet is automatically 
captured by the definition of the capillary number, the 
adsorped nanoparticles cause deviations from the previ- 
ously observed behaviour. At low capillary number and 
low particle coverage, no differences are apparent and 
Taylor's law is reproduced (cf. Fig. [T4|). When the cov- 
erage fraction grows beyond X > 0.40 the deformations 
in this regime increase with increasing X and constant 
capillary number. As it was the case for the system with 
surfactant, the Reynolds number scales linearly with cap- 
illary number. However, since the nanoparticles do not 
affect the surface tension, all curves have the same slope 
(cf. inset of Fig. [Hand Eq. (@H])). Tnis implies that the 
increased deformation in the case of added nanoparticles 
is not caused by changes in inertia of the fluids. On the 
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FIG. 14. Main plot: deformation parameter D as a function 
of the effective capillary number Ca eff for various degrees of 
droplet interface particle coverage fraction x- In an these sim- 
ulations gbr — 0.10, r p = 5.0, m v — 524, and P = 90°. The 
system size is n x = n y = 256, n z = 512. The capillary num- 
ber is varied by changing the shear rate. Lines are added to 
guide the eye and clarify that the effect of adsorped particles 
is very weak for low but that the effect becomes noticeable 
at x > 0.4, where the deformation increases with x a ^ con- 
stant capillary number. In all cases, however, Taylor's law is 
reproduced for small Ca (dashed line). Note that the squares 
in this plot correspond to the triangles in Fig. Inset: as is 
the case with surfactant, the Reynolds number scales linearly 
with the capillary number. Since the nanoparticles do not 
affect the surface tension, all curves have the same slope (cf. 
Eq. ([481Q. 

other hand, the inertia of the particles themselves plays 
a decisive role here. We have investigated the depen- 
dence of the droplet deformation on the size and mass 
of the particles. Particle radii have been varied between 
4.0 < r p < 10.0 and at Ca eff = 0.1 this has led to only 
a small change in D. Yet, changing the mass of the par- 
ticles directly has a substantial effect. We have varied 
the mass of the particles over two orders of magnitude, 
as shown in Fig. [T5j x = 0.55 and Ca eff = 0.1 are kept 
constant and we have rescaled the mass scale with the 
reference mass: m* = m p /524. The particles are ac- 
celerated as long as they are on the part of the droplet 
interface that experiences a shear flow at least partially 
parallel to the particle movement. Eventually, particles 
have to "round the corner" and are forced to move per- 
pendicular to or even antiparallel to the shear flow. The 
increased inertia of heavier particles makes it more diffi- 
cult to change the movement of these particles, leading 
to a situation where the droplet interface is in fact ini- 
tially dragged farther away in the direction of the shear 
flow instead. This process is balanced by the surface 
tension as the surface area increases. This then explains 
the increase of deformation with increasing particle mass. 
As our deformation is increased substantially, the system 



0.45 




0.20 1 , , , ^ 

40 80 120 160 

m* p 



FIG. 15. The deformation parameter D is shown as a func- 
tion of the rescaled mass of the particles m* = m v /m^ where 
rrip = 524 is defined by setting the density of the particle to 1. 
The particles have a radius r p = 5.0, their coverage fraction is 
X = 0.55 and the capillary number is Ca eff = 0.1. Snapshots 
of the droplets are included, showcasing the deformations of 
the droplet. The inertia of the heavier particles causes addi- 
tional deformation as they drag the droplet interface in the 
direction of the shear flow. 

size limits the deformation we can induce. Therefore, the 
values presented here are underpredictions of the actual 
effect of increased mass at high deformations, and might 
indeed hide a breakup event. 

The effect of particles on the inclination angle of the 
droplet is quantified in Fig. [16l We now return to using 
particles with mass m p = 524. At low capillary num- 
ber the disturbance caused to the droplet shape by the 
particles makes the inclination hard to measure. Due to 
the higher Reynolds numbers in these simulations when 
compared to the system with added surfactant, the in- 
clination angle surpasses the 45° mark in all cases, even 
in the case without particles at all |57j . As in the study 
of deformation, the effect of a small number of particles 
is relatively minor, but for x > 0-4, the inclination angle 
decreases sharply, as the droplet becomes more elongated 
and aligned with the shear flow. Increasing the particle 
mass also lowers the inclination angle, for the reasons 
described above, as can be observed in the droplet snap- 
shots in Fig. [15j 



D. Droplet breakup 

When the capillary number is increased beyond the 
values shown in this work so far, we first proceed into a 
regime of extreme droplet deformation, where ellipsoidal 
approximations of the droplet shape no longer hold. This 
is followed by a regime of droplet breakup, where the sur- 
face tension cannot keep the droplet together and two or 
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FIG. 16. Inclination angle 0d as a function of the effective 
capillary number Ca eff for various degrees of droplet interface 
particle coverage fraction In all these simulations gb r = 
0.10, r p = 5.0, m p — 524 and P = 90°. The system size is 
n x = rty = 256, n z = 512. The capillary number is varied by 
changing the shear rate. Due to the higher Reynolds numbers 
in these simulations when compared to the previous system, 
the inclination angle surpasses the 45° mark (dashed line) in 
all cases, even in the case without particles. As in the study 
of deformation, the effect of a small number of particles seems 
to be relatively minor, but for x > 0.4, the inclination angle 
decreases sharply. 



more smaller droplets form. Their increased relative sur- 
face area and smaller volume render them more stable 
against new deformations or breakup events. A series of 
snapshots of this process is shown in Fig. [T3 First, a 
droplet without particles is shown in its steady state (a), 
strongly deformed at an applied shear velocity u s = 0.07, 
but not breaking up. At the same applied shear veloc- 
ity, a particle-covered droplet evolving in time is shown. 
First, deformations take place within the ellipsoidal ap- 
proximation (b & c). As the droplet is deformed even 
more, a definite neck is observed (d & e). When this neck 
pinches off, two droplets are formed. In the highly de- 
formed state just before breakup, the particles are mostly 
found near the centre of the ^-direction, on the parts of 
the interface with highest curvature (this is an extreme 
example of the distributions described in section IIVB[) . 
This means that just after the breakup, even though the 
new droplets are not very strongly deformed, there is a 
large anisotropy in the distribution of the particles, that 
is, one side of each droplet is mostly vacant (f). Af- 
ter more relaxation, however, the particles redistribute 
themselves over the interface in a similar fashion as before 
(g). Analyzing this behaviour in detail remains outside 
the scope of this work. We do remark that introducing 
adsorped particles decreases the resilience of the droplet 
against breakup, effectively lowering the critical capillary 
number at which breakup occurs, which can be viewed 
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FIG. 17. Example of the breakup of a droplet when sub- 
jected to shear flow. The systems shown are identical (n x = 
n y = 256, n z = 512, g br = 0.10, u s = 0.07, R d nit = 76.8, 
Ca eff = 0.15) apart from the introduction of neutrally wetting 
particles of radius r p = 5.0: a) % = 0.0, b-g) x — 0.55. In b- 
g) snapshots of the particle-covered droplet at different times 
are shown. The system without particles has reached a steady 
state at t — 100000 (a). At the same applied shear velocity, 
the droplet with the particles breaks up into two droplets of 
similar size. At b) t = 50000, c) t = 70000 d) t = 80000 and e) 
t — 90000 the droplet still holds together, but the deformation 
is extreme, departing from the ellipsoidal approximation and 
displaying a clear pinch-off. At f) t = 100000 the droplet has 
broken up into two similar-sized droplets, with the particles 
still distributed much as they were on the original droplet. 
After some relaxation the particles have redistributed them- 
selves on the new interfaces at g) t = 150000. 



as an extension of the increased deformations. 



V. CONCLUSION 

In this work we have applied our implementation of 
the lattice Boltzmann method, extended to deal with 
multiple fluid components, surfactants and hard-sphere 
nanoparticles to study physical phenomena related to a 
droplet in shear flow. Surface tensions in a binary system 
can be mapped to the choice of interaction strength be- 
tween the fluid components and can be further adjusted 
by addition of a surfactant species. In this way, the sur- 
face tension can be varied by an order of magitude within 
the stable parameter region. The addition of spherical, 
neutrally wetting particles to the droplet interface does 
not affect the surface tension, owing to the fact that these 
only change interfacial free energy by removing part of 
the energetically unfavourable fluid- fluid interface. 

When a droplet is subjected to simple shear flow, one of 
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the characterizations of the system is the capillary num- 
ber, relating the magnitude of the viscous forces to the 
magnitude of the surface tension. We have found that 
a measured effective shear rate better characterizes the 
system than the imposed shear rate, owing to the distor- 
tion in the velocity fields created by the presence of the 
droplet. 

We have recovered Taylor's law for small deformations 
of a binary droplet, obtaining linear behaviour with the 
analytically predicted slope. For higher capillary num- 
ber, the deformation increases more strongly than this 
linear relation. The surfactant model also conforms to 
this law: when surfactant is introduced into the sys- 
tem the capillary number is changed through the induced 
change in surface tension. Therefore, the same curve as 
found for the binary system is recovered. 

The effect of the addition of nanoparticles adsorped 
to the droplet interface on the deformation properties 
of the droplet has been studied. The particles are not 
homogeneously distributed over the droplet surface, but 
form more densely packed patches in areas with low shear 
velocities and high curvature. This pattern is in a dy- 
namic equilibrium, and particles rotate over the droplet 
interface. Their rotational frequency increases with cap- 
illary number and decreases with distance from the cen- 



tre of the system. For low capillary number or low cov- 
erage of the interface the effect of these nanoparticles 
is negligible. However, in the regime of high capillary 
number and high coverage (~ 50% in the undeformed 
state) the presence of particles induces a larger defor- 
mation at constant capillary number and a decrease in 
inclination angle. This is caused by the inertia of the 
massive particles. Finally, adsorped particles make the 
droplets break up more easily, lowering the critical cap- 
illary number at which breakup occurs. Emulsions con- 
sisting of such particle-covered droplets are expected to 
exhibit shear- thinning behaviour, as the increased defor- 
mation at higher shear rates lowers the apparent viscosity 
of such a complex fluid. 
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